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Abstract 


The authors have been developing an empirical mathematical model to predict the dynamic behaviour of a polymer electrolyte membrane fuel 
cell (PEMFC) stack. Today there is a great number of models, describing steady-state behaviour of fuel cells by estimating the equilibrium voltage 
for a certain set of operating parameters, but models capable of predicting the transient process between two steady-state points are rare. However, 
in automotive applications round about 80% of operating situations are dynamic. To improve the reliability of fuel cell systems by model-based 
control for real-time simulation dynamic fuel cell stack model is needed. Physical motivated models, described by differential equations, usually 
are complex and need a lot of computing time. To meet the real-time capability the focus is set on empirical models. Fuel cells are highly nonlinear 
systems, so often used auto-regressive (AR), output-error (OE) or Box-Jenkins (BJ) models do not accomplish satisfying accuracy. Best results are 
achieved by splitting the behaviour into a nonlinear static and a linear dynamic subsystem, a so-called Uryson-Model. For system identification 
and model validation load steps with different amplitudes are applied to the fuel cell stack at various operation points and the voltage response 
is recorded. The presented model is implemented in MATLAB environment and has a computing time of less than 1 ms per step on a standard 
desktop computer with a 2.8 MHz CPU and 504 MB RAM. Lab tests are carried out at DaimlerChrysler R&D Centre with DaimlerChrysler PEMFC 


hardware and a good agreement is found between model simulations and lab tests. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Vehicle operation is likely to be the most demanding fuel cell 
application. Whereas, the internal combustion engine had more 
than 100 years to meet today’s requirements, fuel cells are forced 
to catch up within a few decades. The main requirements are high 
volumetric and gravimetric power density, reliable operation, 
long lifetimes, low costs and an excellent dynamic behaviour, 
which in total ties a heavy package for fuel cell research. 

Due to their good dynamic behaviour, high power density, 
compact design, non-corrosive electrolyte, simple manufactur- 
ing and a low operation temperature PEMFC are predestined 
for vehicle applications. In automotive applications round about 
80% of operating situations are dynamic. For further develop- 
ment and optimisation of the fuel cell-based power train, models 
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simulating the dynamic behaviour of the system, are of major 
importance. 

In transport applications another model requirement appears. 
Beside describing the dynamics, fuel cell models implemented 
in vehicle controller units must be able to do precise simulations 
in real-time. 

From the theoretical point of view, the voltage of a PEMFC 
stack is a nonlinear multiple input single output (MISO) system. 
The most important variable inputs are the drawn current, fuel 
cell stack temperature, pressure, gas composition and humidity 
at cathode and anode side. For this kind of MISO system, yet 
there is no standardized procedure neither to find a matching 
model structure nor to choose suitable types of models. 

To meet the requirements the focus is set on black-box 
models because physical and physically motivated models like 
Refs. [1] or [2] need a lot of computing time (1h per step 
at Intel Xeon 2.0GHz [1], (100-fold real-time [2]). Unlike 
physical approaches another advantage of black-box modelling 
is the capability to identify system parameters with measured 
data. Using system identification there is no need to know 
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specialized material properties like, e.g. pore size distribution 
of the gas distribution layer or effective active catalyst area. 

Applying a systematical process to system identification 
seems reasonable. The process described in Ref. [3] can be seen 
in Fig. 1. 

Starting with some expert knowledge on the desired system, 
first a sufficient number of experiments have to be designed. 
Mainly, there are statistical and conventional approaches for 
designing experiments. The pros and cons can be read in Ref. [4]. 
Common in most design methods the investigation boundaries 
have to be specified. Some are based on system physics, others 
caused by hardware limiting or safety requirements, respec- 
tively. After this step the desired data has to be examined. The 
quality of the received data will strongly influence the model 
goodness, so the experimental work has to be done conscien- 
tiously. In most cases a post-processing like excepting outliers or 
some signal filtering of the raw data needs to be done for further 
analysis. If the resulting data does not satisfy all quality goals 


no 


satisfying data quality 


Examine the properties of 
the model obtained 


yes 


no 


model good enough 


yes 


model ready for application 


Fig. 1. System identification procedure [3]. 


some details of the experimental design should be improved 
and the following steps must be redone. With high quality data 
a matching model structure has to be selected. One key to suc- 
ceed in system identification based on modelling is to choose an 
appropriate model type. In the following step the model param- 
eters have to be identified. The least squares method is the most 
used parameter identification technique [20]. After the identifi- 
cation step, simulations with the developed model of the desired 
system behaviour can be done and compared to the measured 
data. If the gap between simulated and measured data is within 
acceptable limits, a validation should be done. For validation it 
is important to use a data record, which was not used for param- 
eter identification [3]. If the model does not agree sufficiently 
with the observed data or is not good enough for the proposed 
application another iteration beginning with the model selection 
step must be applied. Can the questions be affirmed, the model 
is ready for application. 

Today modern computing technologies offer various sup- 
ports to handle the model selection, parameter identification as 
well as simulations with the developed model. In these work 
packages the presented work was supported by the widespread 
mathematical software package MATLAB®, provided by The 
MathWorks Inc. 


2. Experimental work 
2.1. Design of experiments 


A conventional design of experiments (DOE) was chosen due 
to low knowledge of the dynamic behaviour of the fuel cell stack. 

The investigation of the static behaviour of the PEMFC stack 
was done at five different current densities (i), three cathode 
excess ratios (A, =ratio of oxygen supplied to oxygen reacted), 
three cathode pressures (P,) and three relative humidity levels 
(th). The anode pressure was adjusted to cathode pressure with 
a constant offset and anode excess ratio was kept at a constant 
level as well as the fuel cell stack temperature. To keep the 
experimental effort acceptable, relative humidity levels at anode 
and cathode were set equal. This leads to 135 measuring points. 

Current dynamics were stimulated by applying positive as 
well as negative current steps each with four different ampli- 
tudes at the same operating point as in case of statically 
characterization. 

Cathode pressure dynamics were also investigated by apply- 
ing pressure steps. In this case transmembrane pressure 
difference has to be limited. This limit depends primarily on 
mechanical properties of the MEA and flow field geometry [5]. 
Cathode pressure steps were applied only with one amplitude 
both in positive as in negative direction related to anode pressure 
level. Therefore, four different steps were investigated. 


2.2. Data examination 


All experimental work was done with DaimlerChrysler lab 
hardware. The investigated stack consists of four cells and 
was operated with pure hydrogen as anode gas and cleaned 
air as cathode gas. Both reactants were supplied via labora- 
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tory infrastructure and flow was controlled by digital mass flow 
controllers. Humidification with deionised water took place in 
heated trickle beds installed in front of the fuel cell stack. 
Pressure adjustment of anode and cathode gas was realised by 
digital controlled valves connected to the air outlet and hydro- 
gen outlet, respectively. Deionised water was used as coolant. 
Temperature was controlled by matching the coolant flow by 
controlling the coolant pump speed. Cathode and anode current 
collector plates were connected to an electronic load via 1 in. 
wires. Parallel to standard data acquisition with a sample rate 
of 1 Hz provided by the test bench fast data acquisition was 
done with a National Instruments PCI DAQ card at a sampling 
rate of 5000 Hz. Test bench operation was controlled by a stan- 
dard desktop computer with National Instruments LabVIEW 
software. 

To ensure comparable experimental setup several important 
parameters were measured. Six thermocouples and pressure 
sensors were used to measure in- and outlet temperatures and 
pressures of cathode, anode and coolant flow. Temperature sen- 
sors were also placed in the humidifiers to guarantee constant 
moisture levels. 

Measuring the 135 statically points first the current density, 
cathode pressure and excess ratio as well as relative humidity 
were set to the desired levels and kept constant for 30 min to 
reach acceptable steady-state conditions. Some quality control 
about reaching steady-state conditions was done by viewing the 
stack voltage time plot. After this 30 min a 2 min mean with a 
resolution of | measurement per second was taken for steady- 
state value. 

Transient voltage behaviour during current variation was 
characterized by standard resolution of 1 Hz over a long time 
period and to get better information about fast current dynamics 
for one second stack voltage was measured with a sampling rate 
of 5000 Hz additionally. 

The response on pressure steps were mapped again with 1 Hz. 

To avoid degradation influence in the recorded data in reg- 
ular time intervals reference VI-characteristics were taken and 
compared to begin of test campaign VI-characteristics. 

Finally, the time elapsed from start to end of simulation was 
measured with a standard watch. 


3. Fuel cell stack model 


System theory provides numerous classes of models each 
with benefits and disadvantages. A brief overview on classes of 
system models is given in Fig. 2. 


3.1. Model structure 


To meet the real-time capability the focus is set on empiri- 
cal models. Fuel cells are highly nonlinear systems [7], so often 
used auto-regressive (AR), output-error (OE) or Box-Jenkins 
(BJ) models accomplish no satisfying accuracy. Best results 
are achieved by splitting the behaviour into a nonlinear static 
and a linear dynamic subsystem, a so-called Uryson-Model 
[6]. 


Fig. 3 shows the developed model structure. 


system models 


linear nonlinear nonlinear 


linear polynoms AR nonlinear DE 
algebraic look-up-tables ARMAX Hammerstein- 
equations ANN OE Model 
fuzzy systems BJ Wiener-Model 
Uryson-Model 
ANN 


fuzz y systems 


Fig. 2. A brief overview on classes of system models (ARMAX, auto-regressive 
moving average model with exogenous inputs; DE, differential equation; ANN, 
artificial neuronal network) [6]. 


3.2. Nonlinear static subsystem 


There are a lot of publications about steady-state behaviour 
of PEMFC [7-14]. 

In this work the focus is set on a dynamic fuel cell stack model 
for real-time simulations. Therefore, an easy to handle and very 
precise element for the nonlinear static subsystem is desired. 
Look-up-tables combine these attributes in the most efficient 
way [21]. 

The empty spaces between the 135 measured points were 
bridged by spline interpolation. In the investigated operating area 
the error between the look-up-table data and the measured stack 
voltages is below 60 mV, so the look-up-table can be considered 
as precise. 

The disadvantage of this method is the relatively great effort. 
If more than four inputs should be varied over a wide range the 
look-up-table method reaches an enormous number of exper- 
iments (e.g. 7 inputs each with 4 levels lead to 47 = 16,384 
experiments!). 
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Fig. 3. Structure of the developed model. 
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Fig. 4. Measured and identified voltage response to a high current density step. 


3.3. Linear dynamic subsystem 


In Ref. [15] some current steps are investigated. The time 
constant of the voltage response is in the range of 50-150 ms. 
Additionally in Ref. [16] a slow process (10-100 s) after current 
steps can be seen. Both observations together suggest a parallel 
structure for current dynamics, one for the slow and another for 
the fast process. 

In Ref. [1] an RC-element models the fast current dynamics in 
a first approximation but is not precise enough. In Ref. [17] there 
can be seen some voltage overshooting at negative current steps 
and some undershooting at positive current steps, respectively. 
A description for this effect is given in Ref. [19]. 

Modelling the transient behaviour in DMFC with transfer 
functions was applied in Ref. [18]. Translated to process models 
a RC-element equals a PT1-element with the transfer function: 


Kı 


Gprı(s) = Geely 


() 
As mentioned just with a PT 1-element the model fit on fast 
current dynamics is not good enough. In some cases the mea- 
sured stack voltage shows a light over- or undershooting (see 
Fig. 4). To improve the model quality parallel a PT2Z-element, 
with the transfer function: 
K2(T,-s+ 1) 


Gpt2z(s) = (s+ 1-5 4+ 7 


was added. 

The slow current dynamics are quite different to the fast 
current dynamics. Here a massive overshoot in case of load 
reduction and an undershoot of stack voltage at positive current 
steps can clearly be seen (see Fig. 7). In the first few seconds 
after a load step a rapid decrease of over- and undershooting can 
be seen. After that initial approach stack voltage comes slowly 
closer to steady-state conditions. Therefore, the slow current 
dynamics were built with three basic elements. Two parallel 
differential elements, so-called DT 1-elements, were connected 
with a static gain element (P-element) in series. 


current voltage 


Fig. 5. Internal structure of the slow current dynamics element. 


The transfer functions of these three elements are: 


en ee 6) 
(T4-s+ 1) 

Gprip(s) = = (4) 
(T5-s+ 1) 

Gp(s) = K3. (5) 


Fig. 5 gives a look inside the slow current dynamics block. 

For the cathode pressure dynamics the same structure as for 
the slow current dynamics was used, because once again over- 
or undershooting of steady-state stack voltage together with a 
faster and a slower approach can be seen. 

In some preliminary investigations gas humidity and cath- 
ode stoichiometry do not show fast dynamic behaviour in the 
investigated operating area. 


4. Simulation results 


The dynamic input signal for each simulation was the mea- 
sured current density or cathode pressure signal, respectively. 
All other model inputs were kept constant. 


4.1. Fast current dynamics 


The shape of voltage response within the first second after a 
load step depends on the operating point, the applied amplitude 
of the step and the direction of the step. 

This reliance prevents the use of fixed model parameters for 
simulating varied conditions. This problem was solved by iden- 
tifying the five current density levels one parameter set for each 
of the eight different load steps in the centre of the operating 
area. 

These 40 records of model parameters were stored in a look- 
up-table. For simulating the fast dynamics of a load steps with 
the amplitude A at a current density B the matching parameter set 
was read from the look-up-table and used for all simulations with 
these AB conditions. If pressure, stoichiometry and humidity 
exactly match the identified case, simulation results fit precisely 
the measured data. If these parameters differ, a certain difference 
between simulation and measured data occurs. Both scenarios 
are shown in Fig. 6. 


4.2. Slow current dynamics 


For simulations a fixed set of model parameters was used for 
the slow current subsystem. The results for four positive and 
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Fig. 6. Comparison of measured and simulated stack voltage in case of fast 
current dynamics (top: all input parameters are the same like in the case used for 
identification, bottom: measured data are at a lower pressure and stoichiometric 
level). 


four negative load steps are shown in Fig. 7. The model shows 
in the centre of the investigated operating area a good agreement 
between simulated and measured fuel cell stack voltage. At the 
edges of the operating area the chosen model parameters do not 
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Fig. 7. Comparison of measured and simulated stack voltage in case of slow 
current dynamics. 
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Fig. 8. Comparison of measured and simulated stack voltage in case of pressure 
dynamics. 


match the measured voltage well. In Fig. 7 this can easily be seen 
at time range between 500 and 700s. Here the highest possible 
current density was applied and the stack voltage drops deeper 
than the model predicts. This failure often can be observed at 
high current densities. 


4.3. Cathode pressure dynamics 


To model the cathode pressure dynamics the same model 
structure as for the slow current dynamics was used. And for 
the pressure model a fixed set of model parameters was used for 
simulations. In Fig. 8 the results of four pressure steps can be 
seen. 

The peak in the pressure line after a step are caused by the 
test bench, especially by the setup of the controller parame- 
ters of the pressure valves. Due to the low sampling rate (1 Hz) 
of this measurement the really occurred height of the pressure 
peaks could be higher than the plotted one. This should be 
taken into account when the quality of the simulation results 
is judged. Once again in the centre of the investigated oper- 
ating area the model fits measured data well apart from the 
very first moments after the step where the goodness is a 
bit lower. At the edges of the investigated area the fit loses 
accuracy. 


4.4. Simulation time 


It took less than 40 s to simulate 100 times a step with a length 
of 1000 ms at a sampling rate of 5000 Hz. This results in a total 
computing time of less than | ms per single step and fulfils the 
requirement. 


5. Conclusion 
With the developed Uryson-Model a modern PEM fuel cell 


stack can be modelled precisely. It is shown that the tran- 
sient behaviour within the first second after current density 
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step can be simulated by a parallel PT1—PT2Z-element. Oper- 
ation point depending model parameters improve the model 
quality in a wide operating area. It is also shown that a 
parallel combination of two DT1-elements in series with a 
P-element both the slow current dynamics as well as the pres- 
sure dynamics can be modelled. Without any optimisations 
in model code or programming work the developed model 
computes dynamic behaviour of a fuel cell stack faster than 
real-time, so the shown process model approach meets this 
requirement. 


6. Future work 


To improve model quality especially at the outer points of 
the operating area in future a mathematical description for the 
dependency of model parameters to the actual operating point 
will be focused. Further dynamic input parameters like humid- 
ity and excess ratio must be added as well as an expansion 
of operating area to reach readiness for automotive applica- 
tions. A reduction of experimental effort should be reached 
by some statistically DOE or an artificial neuronal network 
approach. 
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